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On the basis of near- infrared imaging observations, we derived visual extinction (Ay) 
distribution toward ten Bok globules through measurements of both the color excess 
(Eh-k) and the stellar density at J, H, and K s (star count). Radial column density 
profiles for each globule were analyzed with the Bonnor-Ebert sphere model. Using 
the data of our ten globules and four globules in the literature, we investigated the 
stability of globules on the basis of £ max , which characterizes the Bonnor-Ebert sphere 
as well as the stability of the equilibrium state against the gravitational collapse. We 
found that more than half of starless globules are located near the critical state (£ m ax = 
6.5 ± 2). Thus, we suggest that a nearly critical Bonnor-Ebert sphere characterizes 
the typical density structure of starless globules. Remaining starless globules show 
clearly unstable states (£ max > 10). Since unstable equilibrium states are not long 
maintained, we expect that these globules are on the way to gravitational collapse or 
that they are stabilized by non-thermal support. It was also found that all the star- 
forming globules show unstable solutions of £ max > 10, which is consistent with the 
fact that they have started gravitational collapse. We investigated the evolution of a 
collapsing gas sphere whose initial condition is a nearly critical Bonnor-Ebert sphere. 
We found that the column density profiles of the collapsing sphere mimic those of 
the static Bonnor-Ebert spheres in unstable equilibrium. The collapsing gas sphere 
resembles marginally unstable Bonnor-Ebert spheres for a long time. We found that 
the frequency distribution of £ max for the observed starless globules is consistent with 
that from model calculations of the collapsing sphere. In addition to the near-infrared 
observations, we carried out radio molecular line observations (C 18 and N2H 4 ") toward 
the same ten globules. We confirmed that most of the globules are dominated by thermal 
support. The line width of each globule was used to estimate the cloud temperature 
including the contribution from turbulence, with which we estimated the distance to 
the globules from the Bonnor-Ebert model fitting. 

Subject headings: ISM: clouds — dust, extinction — ISM: globules — stars: formation 
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1. Introduction 

Probing the physical evolution of dense cores toward the onset of star formation is of great 
importance in understanding the nature of the star formation process. The star formation process 
as well as the characteristics of forming stars are probably mostly determined by the properties of 
the nursing cores, such as size, mass, temperature, turbulence, and density structure, right before 
the gravitational collapse. In order to clarify the initial condition of star formation, it is important 
to provide a number of dense core samples with well defined physical properties. 

Bok globules (i.e., isolated molecular cloud cores) are ideal sites for studying isolated star 
formation processes because of their simple shape and isolated geometry from neighboring clouds. 
These characteristics of globules enable us to model their internal structure (e.g., density). For 
these reasons, detailed observations of some well-known globules have been done. Alves et al. 
(2001) reported a high-resolution dust-extinction (Ay) study of the starless globule Barnard 68, 
based on near-infrared observations. They presented a radial column density profile of the globule 
which is well fitted by the Bonnor-Ebert sphere model (Bonnor 1956; Ebert 1955). Because the 
density structure is well determined, molecular abundances in Barnard 68 have been studied by 
comparing dust extinction with molecular distribution (Bergin et al. 2002; Hotzel et al. 2002a, 
2002b; Di Francesco et al. 2002; Lai et al. 2003). 

An unbiased survey of both starless and star-forming globules is important to provide essen- 
tial information on the evolution of density, velocity, and chemical structure of globules as well as 
their relative lifetime. Statistical studies of globules (dense cores) have been conducted based on 
the observations of thermal dust emission (e.g., Visser et al. 2001, 2002), dust extinction in the 
mid-infrared (Bacmann et al. 2000), molecular line emission (e.g., Launhardt et al. 1998; Caselli et 
al. 2002), and a combination of these (e.g., Tafalla et al. 2002; dust and molecular line emissions). 
However, no statistical studies of the density structure of globules based on the near-infrared ex- 
tinction have been carried out. Recently, the dust temperature gradient in cold starless cores (e.g., 
Evans et al. 2001; Ward-Thompson, Andre, & Kirk 2002) has been pointed out, and this may 
systematically affect the column density distribution derived by assuming a constant dust tem- 
perature. The measurements of extinction in the mid-infrared surface photometry (e.g., Bacmann 
et al. 2000) can also be affected by dust temperature distributions, because both absorption and 
emission from the line-of-sight medium should be taken into account at these wavelengths. As a 
column density tracer, near-infrared extinction measurement is one of the most straightforward 
methods with no ambiguity caused by dust temperature distribution. Near-infrared extinction 
measurements can trace dust column densities at low-density outer regions of the globules, where 
(sub)millimeter dust emission is hard to detect. 

The aim of our near-infrared (J, H, and K s ) imaging survey of Bok globules is to reveal the 
evolution of radial density structures of globules from the starless to protostellar phases as well as to 
derive their fundamental physical parameters (e.g., size, temperature, mass, and external pressure) 
on the basis of the measurements of dust extinction. We investigated the radial Ay distribution 
for each globule using the Bonnor-Ebert sphere model. The solutions of the Bonnor-Ebert sphere 
can be characterized by the dimensionless radial parameter £ max , which describes the stability of 
the gas sphere against gravitational collapse, and the critical state is achieved at £ m ax=6.5. By 
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combining our results for ten Bok globules with previous reports on four sources, we investigate the 
difference in physical properties, in particular the £ max value, between the globules in the starless 
phase and those in the star-forming (protostellar or gravitational collase) phase, which is likely 
to be related to the difference in evolutionary stages of the globules. In order to investigate the 
velocity structure of globules, we observed the same globules using the 45 m radio telescope of 
Nobeyama Radio Observatory with the C 18 (J = 1 — > 0) and/or N2H + (J = 1 — ► 0) molecular 
lines. 

The procedures of our observations and data reduction are described in § 2. The distribution 
of dust extinction (Ay) for each globule is obtained in order to derive the radial Ay profile and to 
fit it by the Bonnor-Ebert sphere model (§ 3 and 4). We investigate the difference in the physical 
properties between the starless and star-forming globules and discuss implications for the evolution 
of globules in § 5. Our conclusions are summarized in § 6. 

2. Observations and Data Reduction 

2.1. Source Selection 

We observed ten Bok globules, CB 87, CB 110, CB 131, CB 134, CB 161, CB 184, CB 188, 
FeSt 1-457, Lynds 495, and Lynds 498. These were selected from a unified catalog of dark clouds 
(Dutra &: Bica 2002; compilation of 21 dark cloud catalogs) with the following criteria: (1) globules 
whose angular sizes are less than ~5', (2) those having simple shapes (circular or elliptical), (3) 
those well isolated from neighboring clouds, (4) those located in the solar neighborhood (distance 
of less than ~500 pc), and (5) those observable from both hemispheres. Our selection criteria of 
(1) to (4) are similar to those used in the CB catalog (Clemens & Barvainis 1988), which lists 248 
small, isolated, optically selected dark clouds. Typical angular size and ellipticity of dark clouds in 
the CB catalog are ~ 3' and ~ 2, respectively (see, Figures 1 and 2 of Clemens & Barvainis 1988). 
In the CB catalog, the fraction of large (e.g., 3> 5') and/or highly elliptical (e.g., ellipticity 3> 2) 
clouds with isolated geometry is small. The properties of the observed globules are summarized in 
Table 1. Column 1 lists the name of the globule, Columns 2 and 3 list the equatorial coordinates 
(J2000) of the center of the target field. Column 4 gives the distance to the globule from the catalog 
(Dutra &: Bica 2002). The distance to the nearby and compact globules is difficult to determine. 
The lack of foreground stars toward globules makes spectro-photometric distance measurements 
impossible, and the distance estimation from LSR velocity (kinematic distance) is not reliable 
for nearby clouds. Thus, the listed distances are only values assumed by associating observed 
globules with neighboring larger molecular clouds whose distance is known. Column 5 represents 
the association of IRAS point sources within the optical boundary of each globule. We classified 
globules into starless and star- forming clouds based on the association of IRAS point sources*. A 



*Recently, a low luminosity protostar was discovered toward the globule Lynds 1014, which is known to accompany 
no IRAS point sources, using the Spitzer Space Telescope (Young et al. 2004). Though Spitzer data are useful to 
clarify whether globules contain protostars or not, Spitzer data is available only for one source, CB 131, among the 
"IRAS-less" globules in Table 1. We briefly checked mid-infrared data of the Spitzer telescope toward CB 131 and 
found that there is no protostar candidates near the center of the globule. Though a source is located near the globule 



- 5 - 



Class I protostar (IRAS 19179+1129) is associated with CB 188 (Launhardt 1996; Yun & Clemens 
1992), whereas no evidence of star formation is shown for the other nine globules. For CB 184, a 
Class II source (IRAS 19116+1623) is located near the globule (Launhardt 1996), but its physical 
association with the globule is not known. Since star-forming activity is not taking place around the 
center of CB 184, we will treat CB 184 as a starless globule in this paper. The Bonnor-Ebert model 
fitting studies of globules based on near-infrared extinction have been reported in four sources in 
the literature: Barnard 68 (Alves et al. 2001), Barnard 335 (Harvey et al. 2001), Coalsack Globule 
II (Lada et al. 2004; Racca et al. 2002), and Lynds 694-2 (Harvey et al. 2003). A Class protostar 
(IRAS 19345+0727) is associated with Barnard 335. The other globules have no IRAS sources, 
while Lynds 694-2 shows strong evidence of gas inward motion (Lee, Myers, &: Tafalla 1999, 2001; 
Lee, Myers, & Plume 2004). Thus, we will classify Lynds 694-2 as a star- forming globule in this 
paper. There are two recent reports on the density structure of Coalsack Globule II based on similar 
analyses. We mainly refer to the observations of Lada et al. (2004) because their near-infrared 
observations are more sensitive than those of Racca et al. (2002). 



2.2. Near-Infrared Data 

We carried out near-infrared (J, H, and K s ) imaging observations of Bok globules using the in- 
frared camera SIRIUS (Simultaneous three-color InfraRed Imager for Unbiased Surveys; Nagayama 
et al. 2003; Nagashima et al. 1999) on the 1.4 m telescope, IRSF (InfraRed Survey Facility), at 
the South African Astronomical Observatory (SAAO). Observations were made over two periods 
from 2002 July 30 to August 12 and from 2003 June 24 to July 14. SIRIUS has three 1024x1024 
pixel HgCdTe infrared detectors (HAWAII array), which enables us to take J, H, and K s images 
simultaneously. IRSF/SIRIUS can image a large FOV of ~7'7x7'7 with a scale of 0"45/pixel, and 
the limiting magnitudes (S/N=10) for a typical 10-minute exposure time reach 18.9, 18.3, and 17.3 
mag at J, H, and K s , respectively. Thus, this survey is much deeper than previous extensive sur- 
veys (e.g., 2MASS, DENIS). For each globule, we performed a 15-minute integration with a set of 
dithered 20 s or 30 s exposures except CB 110 and FeSt 1-457, whose integration times are 45 and 
60 minutes, respectively. The dithered spatial interval was 30" or 40". Since the apparent diameter 
of each target globule is less than 5', we covered its entire extent at one time. The typical seeing 
during the observations was 0"9— 1"4 (FWHM), corresponding to 2—3 pixels. To correct detector 
flatness and pixel-to-pixel variation, we obtained twilight flat and dark frames at the beginning and 
end of the nights. We observed some sets of near-infrared standard stars listed by Persson et al. 
(1998) during the observations for photometric calibration. 

We processed the observed data in a standard way of near- infrared image reduction (i.e., dark 
subtraction, flat-field correction, and median sky subtraction) with the IRAF^ (Image Reduction 
and Analysis Facility) and the IDL (Interactive Data Language, by Research Systems, Inc.) rou- 
tines. After subtraction of averaged dark frame, each object frame was divided by normalized flat 



on the 24 fira image, it is not known whether the source is a protostar associated with CB 131. 

'IRAF is distributed by the National Optical Astronomy Observatories, which are operated by the Association of 
Universities for Research in Astronomy, Inc., under cooperative agreement with the National Science Foundation. 
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frames. For correcting thermal emission pattern and the fringe pattern from the atomospheric OH 
emission, we constructed a median sky frame to be subtracted from the object frame. Since all 
the observed targets are located in the very crowded field of the Milky Way, we did not adopt 
the self-sky, with which the residual signal of (bright) stars cannot be excluded on the resulting 
self-sky frame. Therefore, we constructed the sky frame for each band by taking the median of 
dithered frames toward the different several (four to eight) objects observed in the same night, 
which resulted in a sky frame with no stellar residuals. This "average" median sky does not re- 
flect the actual sky for each object because the median was taken across several object frames. 
We compared the "average" median sky frame with the self-sky frame for each object in order to 
correct their different sky levels, and the final (level-matched) median sky frames were made for 
each object. Though the thermal emission pattern is adequately corrected in the final median sky 
frames, the OH fringe pattern is time-averaged. We confirmed that the effect of the time-variation 
of the OH fringe pattern is within acceptable levels by checking the sky noise level of the final 
combined object frames. After the sky subtraction and frame-registration, we combined all the 
dithered images for each object with a 3a clipping. We show three-color (JHK S ) composite images 
of globules in Figure 1. 

For the source detection and photometry, we applied the IRAF/DAOPHOT package (Stetson 
1987) to the images. We measured the instrumental magnitudes of the detected stars having a 
peak intensity greater than 4.0 a noise level above the sky background by fitting a point-spread 
function (PSF) and subtracting it from the image. We iteratively performed stellar detections and 
psf-fitting photometries to the source-subtracted image in order to collect overlooked sources in 
the previous procedure. We manually checked the results of the photometries and removed some 
false detections by visual inspection. We summarize the results of stellar detection in Table 2. For 
calculating the plate solution, we matched the pixel coordinates of a number of detected stars with 
the celestial coordinates of their counterparts in the 2MASS point source catalog (Cutri et al. 2003) 
and obtained the transformation coefficient using the IRAF/IMCOORDS package, which resulted 
in an acceptable positional error of less than 0''5 (rms). All the stars detected in our observations 
were cataloged separately with respect to the filters (JHK S ). These stars of single-band detections 
were used to generate an extinction map at each band with a star count method. We applied the 
following equations for the color conversion from the IRSF to the standard CTIO/CIT photometric 
system: 



These were derived from measurements of red standard stars (see, Nakajima et al. 2004 for details). 
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(1) 

(2) 
(3) 



2.3. 



Molecular Line Data 



We have carried out molecular line observations toward ten globules using the 45 m radio 
telescope of Nobeyama Radio Observatory over four observing periods (2003 March 10-19, 2003 
April 8-10, 2004 January 9-12, and 2004 February 3-10). The observed objects are the same as those 
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observed by IRSF/SIRIUS. The receiver employed was the 25-element focal-plane array receiver 
BEARS consisting of double-sideband SIS mixers (Sunada et al. 2000). All the target globules in 
Table 1 are mapped in position-switching mode in C 18 (J = 1 — ► 0) at 109.782173 GHz (JPL 
catalog; see, http://spec.jpl.nasa.gov/) and/or N2H + (J = 1 — > 0) at 93.1737767 GHz (Casclli, 
Myers & Thaddeus 1995) on a 13" 7 grid spacing. The half-power beamwidth for each element 
beam was estimated to be ~17"8 at 93 GHz and ~15"1 at 109 GHz. The map coverage (~4'x4') 
is sufficiently large compared with the extinction feature of globules. As a back end, we used a 
digital autocorrelator with a 37.8 kHz frequency resolution (~0.12 kms -1 at 93 GHz, and ~0.10 
kms -1 at 109 GHz). The intensity was calibrated using the standard chopper-wheel method, and 
is expressed in terms of the corrected antenna temperature T A . To correct for the sideband ratio, 
we calibrated T* A using the intensity obtained with the single-sideband receiver S100. The scaling 
factor from T* A of BEARS to that of S100 for N 2 H+ (J = 1 — > 0) was measured by ourselves, and 
that for C 18 (J = 1 — > 0) was provided by the observatory. A typical value of the scaling factor 
for each BEARS element was ~ 1.5. We estimated the main beam efficiency of S100 at 93 and 109 
GHz to be 0.515 and 0.442, respectively, by interpolating values at 86, 100, and 115 GHz provided 
by the observatory. The pointing accuracy was typically better than 7" by observing SiO maser 
sources at 43 GHz every 1-1.5 hours during the observations. 



3. Derivation of Ay and Molecular Line Parameters 

Ay measurement is one of the most straightforward ways to probe the density structure of 
dark clouds. There are two standard ways to do extinction mapping; one is the star count method 
(e.g., Dickman 1978; Cambresy 1999; Dobashi et al. 2005) based on the measurements of the stellar 
density, and the other is the near-infrared color excess (NICE) method (Lada et al. 1994) based on 
measurements of the stellar color excess (e.g., Eh-k)- By considering the stellar density in H and 
K s at the densest region of globules, we decided which method is better for each target globule. We 
analyzed CB 131, CB 134 and CB 188 with the star count method and analyzed the other globules 
with the NICE method. 



3.1. Star Count 

On the basis of the stars detected in the J, H, and K s bands, we derived extinction maps (Aj, 
Ah, and Ak s ) for CB 131, CB 134, and CB 188 using the star count method (see, e.g., Kandori et 
al. 2003). A logarithmic cumulative stellar density logiV measured at the magnitude m^o can be 
converted to extinction A\ as 

A x (Aa, AS, m A)0 ) = m A)0 - f' 1 [log N (Aa, A<5, m x ,o)] (4) 

where Aa and A5 are the offset positions from the equatorial coordinate of image center, and 
f^ 1 is the inverse function of / = log 7V re f {m\), which is the logarithmic cumulative luminosity 
function (i.e., Wolf diagram; Wolf 1923) constructed in the extinction- free reference field (see Fig. 
2). Though / is often assumed to be linear as / = a + bm\, we fitted it with a 4th-order polynomial 
because the function shows a slightly curved shape. 



-8- 



To measure stellar density distribution, log N(Aa, AS, m^o), we arranged circular cells (typi- 
cally ~30" in diameter) separated by 9" intervals on each observed field (~7'x7' area), and counted 
the number of stars falling in the cells up to the threshold magnitude m^o (see Table 2 for the 
definition) which is well above the limiting magnitudes. The diameter of circular cells depends on 
the stellar density in each cloud field. We smoothed the stellar density map with a Gaussian filter 
(FWHM=13''5) to reduce noise. In order to derive Aj, Ah, and Ak s , we numerically calculated 
f^ 1 in equation (4) after substituting the actually measured stellar density log N(Aa, A6, m\p) 
into In the resulting extinction maps, Aj tends to be saturated at the densest region of the 
cloud, but sufficiently reveals the diffuse extinction feature around the cloud, and Ak s has the 
opposite tendency. In order to relieve saturation, we made a composite Ay map by combining the 
Aj, Ah, and Ak s data with the following equation: 

A v = aR KsV A Ks + (1 - a)R HV A H for Ah > ^4^ ax /2 (5) 
Ay = (3R HV A H + (l-(3)R JV Aj for Ah < A^ ax /2 (6) 

where a and (3 are the weighting parameters defined as An/A H ax and An/(A H ax /2) where A H &X 
denotes the maximum value of Ah in each analyzing field. The R\y denotes the standard conversion 
factor from A\ to Ay in Rieke & Lebofsky (1985), which is 3.55 and 5.71 for J and H, respectively. 
For K s , we assumed Ay = QAAAk 3 on the basis of equation (3) and the conversion factors in Rieke 
& Lebofsky (1985). These conversion factors depend on the assumption of dust optical properties, 
in particular Ry. We discuss systematic uncertainties in the derivation of Ay in the Appendix. The 
resulting Ay maps for CB 131, CB 134, and CB 188 are shown in Figure 2. Typical uncertainty of 
Ay in the reference field (Ay ~ mag) is estimated to be ~0.9 mag. 



3.2. Near-Infrared Color Excess 

On the basis of the stars detected both in the H, and K s bands, we derived extinction (Ay) 
maps for CB 87, CB 110, CB 161, CB 184, FeSt 1-457, Lynds 495, and Lynds 498 using the NICE 
method. If we assume that the population of stars toward the observed field is invariable, the mean 
stellar color in the reference field < H — K > re f can be used as the zero point of color excess. The 
color excess distribution toward globules can be derived by subtracting < H — K > rc f from the 
observed H — K of stars. Since the extent of each analyzed field is relatively small (~ 7' x 7'), our 
assumption on invariant stellar population seems plausible. To measure color excess distribution, 
we arranged circular cells (~30" diameter) spaced at 9" (the same alignment as used in the star 
count analysis) on each observed field. We calculated the mean stellar color in each cell and derived 
their color excess as 

(A a ,AS)-{H-K) ref , (7) 

where (H — K)i is the color index of the i-th star in a cell, and N is the number of stars falling in 
a cell. (H — K) Te { is the mean color of stars in the reference field (see, Fig. 2). We converted the 
H—K color excess to Ay using the reddening law of Rieke & Lebofsky (1985) as Ay = 15.9 x Eh-k- 
We then smoothed the resulting Ay map in the same manner as used in the star count analysis. 



E H -k (Aa, AS) = 



N 

E 

.i=i 



N 
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The resulting Ay maps are shown in Figure 2. Typical uncertainty of Ay in the reference field 
(Ay ~ mag) is estimated to be ~0.6 mag. 



We derived the observed molecular line parameters, LSR velocity, line width, and optical depth, 
for each globule. We observed N 2 H+ (J = 1 — > 0) for FeSt 1-457 and C 18 (J = 1 — ► 0) for the 
other globules. We show the spectrum observed at the center of each globule in Figure 3. For the 
N 2 H+ spectrum fitting we took into account their hyperfine structure and line optical depth effect 
(using the same procedure as described in Tatematsu et al. 2004). The intrinsic relative intensity 
and rest frequency of hyperfine components are taken from Tine et al. (2000) and Caselli et al. 
(1995), respectively. 

The observed line width AV^bs of the globule ranges from ~ 0.2 to ~ 0.5 km s _1 (Table 3). The 
kinetic temperature T k of low mass dense cores is known to be ~ 10 K (e.g., Benson & Myers 1989). 
Thus, it is natural to attribute the excess in the observed line width with respect to the thermal 
line width (0.13 km s~ x for C 18 and N2H + ) to the contribution from non-thermal turbulence in 
globules. We estimated non-thermal line width AVnt of each globule by assuming that T k is equal 
to 10 K. AVnt can be obtained from AVJ| T = AV^ S — 81n2fcT k /m b s , where m b s is the mass of 
the observed molecule. The thermal line width AVt = (81n2A;T k /m) 1 / 2 for T k = 10 K is 0.443 km 
s -1 , where m is the mean molecular weight (2.33 amu). Derived line parameters are listed in Table 
3. We found that thermal line width for the mean molecular weight is larger than non-thermal 
line width for most of the observed globules. It is likely that observed globules are dominated by 
thermal support with smaller contribution from the non-thermal turbulence. The thermal and non- 
thermal components of the line width can be characterized by a temperature as T k = AV^m/fc81n2 
and Tnt = /^■V^ T m/k8\n2, respectively. We derived T e g = T k + Tnt for each globule in Table 3 
as a cloud-supporting temperature which includes non-thermal contribution. In order to estimate 
the distance to the globules, we will use T e g in combination with the Bonnor-Ebert model fitting 
of globules in §4.2.2. In this paper, we only use the information on molecular line width. Detailed 
studies such as comparison between dust extinction and molecular distribution will be reported in 
a subsequent paper. 



In order to model the internal density structure of Bok globules, we used the Bonnor-Ebert 
sphere model (Bonnor 1956; Ebert 1955), which is a modified Lane-Emden equation describing a 
pressure-confined, self-gravitating isothermal gas sphere in hydrostatic equilibrium. The density 
structure of a Bonnor-Ebert sphere can be obtained by solving the following equation: 



3.3. 



Molecular Line Parameters 



4. Density Structure 



4.1. Bonnor-Ebert Sphere Model 




(8) 



-10- 



where 4> is the logarithmic density contrast, </>(£) = — ln(p/p c ); p and p c are the volume density and 
central volume density, respectively. £ is the dimensionless radial parameter, £ = (r/C s )v / 4vrGp c ; 
r, G, and C s are the radius, gravitational constant, and isothermal sound speed, C s = (kT/m) 1 ^ 2 , 
where k, T, and m are the Boltzmann constant, kinetic temperature, and mean mass of the molecule 
(2.33 amu), respectively. Under the standard boundary conditions, </>(0) = (i.e., p = p c at r = 0) 
and d<j)(0)/d£ = (i.e., dp/dr = at r = 0), equation (8) can be solved with numerical integration. 
If the gas sphere is confined by external pressure P e xt at the core boundary radius R, the solutions 
of the Bonnor-Ebert sphere can be characterized by the dimensionless radius £ max = £(r = R) as 

R 



U = yrV^Gpc. (9) 

£max is a stability measure of the gas sphere against the gravitational collapse (Bonnor 1956; Ebert 
1955). Critical state is achieved at £ max = 6.5 corresponding to the density contrast of p c / Pr = 14. 
For a solution with £ max > 6.5, the equilibrium state is unstable to the gravitational collapse. The 
Bonnor-Ebert sphere mass M and the external pressure P ext are calculated as 

^ max 

P cx t = C s Ve- 0(5max) (11) 



where a = C s / (AirGp, 



4.2. Bonnor-Ebert Model Fitting of Bok Globules 

4-2.1. Derivation of Column Density Profile 

For a comparison of observational data (i.e., Ay) with a theoretical Bonnor-Ebert density 
distribution, radial Ay profiles for each globule were constructed. We set annuli spaced at 9" 
around the center of the core that is the Ay intensity peak (see Table 4), and calculated averaged 
Ay values at each annulus (see Figure 4). For apparently elliptical globules (CB 87, CB 161, and 
CB 184) we performed ellipse fitting to the Ay distribution, and "equivalent" radius, r = Vab, 
was substituted for radius where a and b are the respective semi-minor and semi-major axes of the 
ellipse for the cores. Fitted ellipticity and position angle for the globules are listed in Table 4. It 
is possible that the observed globules with circular projected shapes are actually pole-on prolate 
ellipsoids. In this case, the shape of the observed Ay profile should be steeper than that for the 
intrinsic shape. Since there is no effective method for estimating how ellongated a globule is along 
the line of sight, we assumed that the line of sight length of globules is the same as the projected 
diameter. 

We manually masked extinction features likely to be unrelated to the globules (e.g., the other 
clouds on the same image, diffuse streaming feature; see Fig. 2), and masked regions were omitted in 
the circularly averaging procedure. Furthermore we recalibrated the relative zero point of extinction 
for the analyzed globules by subtracting background extinction (i.e., Ay offset value). The Ay 
offsets have been derived from averaging Ay within the annulus of 150" to 200" radius, which is 
sufficiently larger than the primary extinction feature of globules. Actually, the subtracted values 
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as Ay offset are small (typically ~ 0.5 mag), and it solely depends on where we set the reference 
fields in the star count or color excess analysis. In order to convert Ay to the H2 column density, 
we used the relationship between the optical color excess Eb-v an d the hydrogen column density, 
N(HI + H 2 ) /Eb-v = 5-8 x 10 21 cm" 2 mag" 1 (Bohlin, Savage, & Drake 1978), where N(Hl + H 2 ) 
= TV (HI) + 2iY(H 2 ). This relationship is based on the UV absorption line measurements toward the 
stars of slightly obscured lines of sight (Eb-v ^ 0.5 mag). By assuming that (1) all hydrogen in the 
cloud is in a molecular form and that (2) the ratio of total-to-selective extinction l?y(= Ay / 'Eb-v) 
is equal to the standard ISM value of 3.1 (e.g., Whittet 1992), we obtained N(B. 2 )/A V = 9.4 x 10 20 
cm -2 mag -1 , and used it to derive the H2 column density. The resulting column density profiles 
are shown in Figure 4. 



4-2.2. Derivation of Model Parameters 

In order to produce Bonnor-Ebert model density profiles, we performed numerical integration 
of equation (8). Since the Bonnor-Ebert sphere is characterized by a single parameter £ max (Eq. 
[9]), best-fit £ max for each globule can be determined by comparing model column density profiles 
with those obtained from observations. We integrated the model density profile along the lines of 
sight on grids to calculate the column density profile. In order to match the resolution between 
the model and the observed column density data, the model column density on the Aa-A5 plane 
was convolved with the beam used in the star count or color excess analysis (cylinder-shaped beam 
convolved with the Gaussian). In the x 2 -procedure, the radial profile was fitted from the minimum 
radius (center) to the radius where the observed Ay value reaches ~ 1 mag (detection limit). In 
some cores (Lynds 495 and Lynds 498), deviation of observed profiles from the Bonnor-Ebert sphere 
is found at diffuse outskirts (Ay < 1 mag), which was excluded in the fitting. The x 2 -fitting results 
are not sensitive to the selection of the maximum radius for fitting part unless there is no significant 
kink in the column density profile. 

In Figure 4, we show the results of the Bonnor-Ebert fit to the radial column density profile 
of the globules. The dots and error bars represent the average iV(H2) values at each annulus of 9" 
intervals and the rms deviation of data points in each annulus, respectively. The solid line denotes 
the radial column density profile of the best-fit Bonnor-Ebert sphere which is convolved with beam 
used in the Ay measurements. The dashed line denotes the Bonnor-Ebert model profile before 
the convolution. In Table 4, we list the derived Bonnor-Ebert parameters for globules including 
previous reports on four sources (Barnard 68, Barnard 335, Coalsack Globule II, and Lynds 694-2). 
The error of each parameter in the table represents la deviation in the x 2 density profile fitting. 
We discuss possible systematic uncertainties in the derivation of £ max in the Appendix. 

Recently, Hotzel et al. (2002b) and Lai et al. (2003) presented detailed discussion on the 
relation between the physical parameters of the Bonnor-Ebert sphere. Substituting the relations of 
R = OrcL oc d and p c oc iV^H^)^ 1 ^ -1 oc d" 1 into equation (9), we can obtain d -1 T =constant (Lai 
et al. 2003). The distance d or temperature T of globules should be accurately determined from 
observations. Since there were no accurate methods for measuring distance to the starless, compact, 
and nearby (d < 500 pc) globules except for the limited cases (e.g., interaction with HII region of 
known distance, existence of foreground stars toward them), the distances listed in Table 1 serve 
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as the first reliable values. If we use the initially assumed distance d[ n i, the necessary temperature 
for the Bonnor-Ebert sphere Tbe can be determined. If we assume that T e g(= + Tnt) from 
our molecular line observations (see §3.3) can approximate effective supporting temperature of the 
cloud^, the corrected distance to the globule, d = (T e ff /Tbe) dmi, can be determined. Finally we 
derived all the physical parameters of best-fit Bonnor-Ebert spheres (e.g., radius, central density, 
and mass) for the above "distance-assumed" and "temperature-assumed" cases as listed in Table 5. 
We note that the "temperature-assumed" distance to the starless globule FeSt 1-457 is estimated 
to be ~ 70 pc, which would make this globule one of the nearest dark clouds. 

5. Discussion: Stability and Evolution of Globules 

In this section, we discuss the internal structure and stability of globules by using the data of 
the ten globules from our observations and four from the literature. 

5.1. Physical Properties of Globules 

Though we have estimated physical properties of globules by assuming the Bonnor-Ebert 
model, there are actually other theoretical models (e.g., Shu 1977; Larson 1969; Penston 1969; 
McLaughlin & Pudritz 1997; Plummer 1911). Bacmann et al. (2000) presented that a Bonnor- 
Ebert-like model (a finite size sphere with inner uniform density region and p oc r~ 2 outer envelope) 
fits well with the column density profile of globules obtained from mid-infrared 7 /im absorption 
measurements, but that a singular-isothermal-sphere (SIS) or singular logotrope sphere cannot 
reproduce the observed profile well. By considering the resolution of our extinction maps (~30") 
and possible geometric error of globules (e.g., deviation from spherical symmetry), it is not possible 
to make meaningful comparisons between the Bonnor-Ebert model and other Bonnor-Ebert-like 
models (see, e.g., Figure 2 in Harvey et al. 2003, which compares the Bonnor-Ebert model with a 
Plummer-like model). Here, we simply assume the Bonnor-Ebert model. 

As described in the previous section, the dimensionless radial parameter £ max determines the 
shape of a Bonnor-Ebert density profile as well as the stability of the equilibrium state against 
the gravitational collapse. We plot globules on the £ max versus density contrast diagram in the 
upper panel of Figure 5. Since £ max has one-to-one correspondence with density contrast, all the 
Bonnor-Ebert spheres are distributed along the solid curved line. We found that more than half 
of starless globules (7 out of 11 sources) are located near the critical state, £ max = 6.5 ± 2. Thus, 
we suggest that a nearly critical Bonnor-Ebert sphere characterizes the typical density structure of 
starless globules. The remaining starless globules show clearly unstable states (<£max 

> 10). When 

we divide starless globules into two groups with respect to the critical line (£ m ax = 6.5), there 
are three stable starless globules and eight unstable starless globules. The majority of starless 



* Theoretically, supporting force in a Bonnor-Ebert sphere is purely thermal pressure, but actual clouds also have 
non-thermal support, e.g., turbulence and magnetic fields, which are widely observed in molecular clouds (see, e.g., 
Larson 2003). We assume that non-thermal support can be included in the model by replacing the isothermal sound 
speed C s = (fcTk/m) 1 ' 2 in equation (9) by the effective sound speed C c s = (fcT ff /m) 1 ' 2 . 
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globules is located in the unstable states (£ max > 6.5) if the uncertainties in £max are not taken 
into account. We also found that all the star-forming globules have larger £ max values (> 10) than 
that of the critical Bonnor-Ebert equilibrium state, which is consistent with the fact that they have 
started gravitational collapse. In the lower panel of Figure 5, we show the fitting error in £ max 
for each globule. Filled gray circles denote the fitting results for some globules with considerable 
features of ambient extinction if the masked region on the Ay map (see Fig 2) is included in the 
derivation of the column density profile. In addition to the £ max diagram in Figure 5, we made 
histograms of logarithmic density contrast for starless and star- forming globules (Figure 6) . Similar 
characteristics as descrived above can also be seen in the histograms. 

To see the relationships between derived physical properties of globules in Tables 4 and 5, we 
made correlation diagrams between density contrast and the other physical parameters of globules 
in Figure 7. We consider that the horizontal axis (density contrast) represents the evolutionary 
states of globules, because globules should evolve toward higher central condensation to form stars. 
We can see from panels (e) and (f) in the figure that temperature T e g and external pressure P cx t 
of globules appear roughly constant regardless of the density contrast. If we assume constant T e ff 
and -Pext; the Bonnor-Ebert model parameters for each density contrast value can be constrained. 
The broken line in each panel is the relationship for the Bonnor-Ebert spheres with a constant T e g 
of 13.2 K and P e xt of 5.3 x 10 4 K cm -3 , which are the median values for the globule samples. The 
value of -P cxt is larger than the previously reported ISM pressure in the solar vicinity of ~ 1.8 x 10 4 
K cm" 3 (McKee 1999). 

The overall trend in the data point distributions in Figure 7 can be approximated by Bonnor- 
Ebert spheres with a constant "effective" temperature and a constant external pressure. The 
relations for the radius and mass (panels [c] and [d]) remain roughly constant against the density 
contrast, suggesting that these quantities are not sensitive to the evolutionary states of globules. 
This result agrees with the N2H + observations of molecular cloud cores in Taurus (Tatematsu et 
al. 2004), in which there is no remarkable difference in core radius and mass between starless cores 
and cores with stars. Since the density at core boundary is derived from P cx t/T c s, the density 
at core center n c is directly obtained for each density contrast value and shows a proportional 
relationship against the density contrast as shown in panel (b). Accordingly, the density contrast 
and Ay show positive correlation. The data points in panel (a) are located slightly lower than 
the broken line because of a beam dilution effect in the Ay measurements with ~ 30" resolution. 
For molecular cloud cores located in regions of higher pressure (e.g., dark cloud complex) than in 
the solar vicinity, the Bonnor-Ebert relationship in panels (a) and (b) should be shifted upward 
corresponding to the larger value of P cxt /T e gf. We can see from panel (a) that the density contrast 
value is simply estimated from the Ay (column density) measurements under the assumption of 
the Bonnor-Ebert sphere with a constant T e g- and P e xt- 



5.2. On the Selection Bias of the Globule Samples 

Here we briefly consider the selection bias in our globule samples. The peak Ay values for our 
observed globules range from ~ 6.5 to 41 mag. We selected globules that were opaque at optical 
wavelengths, i.e., with few stars in the background of the cloud on the Digitized Sky Survey image. 
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These sources are expected to have Ay greater than at least 5 mag, and expected to be sites of 
future star formation, i.e., so-called "dense cores". The lower limit of the peak Ay for the globules 
may be determined by the criterion in our target selection. Thus, it is possible that the diffuse cores 
of low column densities are not included in our globule sample as in the case of previous studies on 
dense cores (e.g., Visser et al. 2002; Lee h Myers 1999; Jijina et al. 1999). We note that the lower 
limit Ay of ~ 5 mag corresponds to a low density contrast of ~ 3 if we assume a constant T e g and 
P C xt used in the previous section. The density contrast value is much less than the critical value of 
14. 



5.3. Implications for the Stability of Globules 

It is interesting that there are starless globules showing larger density contrast than that of the 
maximum value (~14) for stable Bonnor-Ebert spheres (see § 5.1), because unstable equilibrium 
states should not last long without any extra stabilizing force such as magnetic pressure and/or 
turbulent pressure (e.g., McKee & Holliman 1999). Since unstable equilibrium states cannot last 
long, the globules of large density contrast should not exactly be the Bonnor-Ebert sphere even 
if they mimic the Bonnor-Ebert density profile. Similar results have also been reported from the 
mid-infrared extinction study of starless dense cores (density contrasts are generally ~10— 80; see 
Table 3 of Bacmann et al. 2000). If we consider a simple case in which the non-thermal pressure 
components have the same radial dependence as the thermal pressure, the theoretical Bonnor-Ebert 
model can be modified by replacing its kinetic temperature term with T e g as described in § 4.2.2. 
In this case, the shape of the Bonnor-Ebert density profile as well as the stability criteria (i.e., the 
value of £ max or maximum density contrast) hold unchanged compared with the purely thermally 
supported case. This modification cannot provide a realistic explanation for the stability of the 
large density contrast globules, suggesting either that they are already collapsing toward higher 
central condensation as discussed in § 5.5, or that the stabilizing mechanisms accounting for their 
large density contrast are more complex. The Bonnor-Ebert assumption for the globules is not 
valid in these cases. We note that our result is in marked contrast to one in which most of the 
clumps observed at 850 fim dust emission in the Orion B molecular cloud complex are in stable 
equilibrium based on the Bonnor-Ebert analysis (Johnstone et al. 2001). 

A non-isentropic multi-pressure polytrope model (McKee & Holliman 1999) may provide suc- 
cessful interpretation allowing larger density contrast than that of the critical Bonnor-Ebert sphere. 
Curry & McKee (2000) developed a composite polytrope model which also allows large density con- 
trast. These two models provide a stable equilibrium solution of arbitrarily large density contrast 
by adjusting the contributions from the non-thermal pressure components, i.e., turbulent and mag- 
netic pressure. Galli et al. (2002) studied the effects of the interstellar radiation field (ISRF) on the 
cloud stability. They reported that a cloud's maximum density contrast increases with increasing 
intensity of the ISRF if the stabilizing effect of external heating of the cloud as well as produced 
temperature gradient are considered. These models are possible alternatives to the Bonnor-Ebert 
sphere providing larger value of maximum density contrast (> 14). The shape of their density 
profiles have characteristics similar to those of the Bonnor-Ebert sphere (see, Fig. 3 and 4 in Galli 
et al. 2002, Fig. 10 in McKee & Holliman 1999), that is a flat inner region surrounded by steeper 



-15- 



envelope decreasing with power-law, in particular p oc r~ a where a ~ 2. It seems difficult to 
distinguish them only from the observed shape of the density profile. In other words, each model 
should provide similar density profiles and corresponding density contrast values for the same ob- 
servational dataset. Detailed observations of velocity and temperature structure of dense cores are 
needed to determine the most realistic model. 

Major differences among the models discussed above are the supporting mechanisms of the 
cloud and the stability criterion, i.e. the value of the critical density contrast. The actual density 
contrast value for globules can be obtained through the Bonnor-Ebert fit (see Table 4), but its 
critical value is not always known (i.e., model dependent). Initial density structure of globules 
right before the collapse can be determined if we constrain the value of its critical density contrast. 
Initial density conditions for star formation will provide a good starting point for the theoretical 
calculations. For example, given initial density profile of a globule, time evolution of mass accretion 
rate can be calculated (e.g. Whitworth & Ward-Thompson, 2001), which is an important parameter 
for constraining the mass of a newly formed star. 

In order to estimate the value of critical density contrast from observations, it should be worth 
investigating the frequency distribution of density contrast for a number of starless globules. Since 
unstable equilibrium states should not be long sustainable, we expect that the detection probability 
for stable globules is much higher than that for unstable ones. Thus, it is likely that the frequency 
distribution of density contrast shows a sharp decrease toward larger density contrasts at a certain 
bin, which should correspond to the density contrast critical for collapse. Though the number of 
our samples (11 sources) is not large, we present the histogram of density contrast for the starless 
globules in Figure 6. The density contrast of the starless globules is mostly populated near the 
critical Bonnor-Ebert sphere and decreases toward the larger density contrast, suggesting that the 
actual critical density contrast of globules is not very different from that for the Bonnor-Ebert 
sphere. A statistical study with a larger number of samples is our future plan. 

5.4. Starless Globules with Stable Bonnor-Ebert Solutions 

In Table 4, three out of fourteen globules are fitted as a stable equilibrium state of the Bonnor- 
Ebert sphere. If they are truly stable, how can stable globules form stars? Stable globules might 
increase their density contrast quasi-statically toward the onset of dynamical collapse with gradually 
increasing external pressure (e.g., Hennebelle et al. 2003) and/or with decreasing internal turbulent 
pressure (Nakano 1998). For the case of rapid increase of external pressure, core collapse triggered 
by a compression wave should be taken into account (e.g., Hennebelle et al. 2003, 2004; Motoyama 
et al. 2004). Our observations of globules, however, show neither systematic decrease in T e g- nor 
increase in P cxt (panels e and f of Fig. 7) toward larger density contrast. It seems that non-thermal 
turbulence is already mostly dissipated even in the stable globules. Though there is only small 
room for the dissipation of turbulence, a slight decrease of turbulence may initiate the collapse of 
stable globules as suggested by Nakano (1998). The decrease of T e g (i.e., line widths) necessary for 
stable globules to start gravitational collapse is small. This is consistent with the result that there 
is no apparent tendensity of systematic decrease in T c g against density contrast (panel e of Fig. 
7). In the star forming regions where non-thermal motions are more prominent than globules, e.g., 
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GMCs, the dissipation of turbulence in dense cores is reported through the virial analysis (e.g., 
OMC-2/3 region in Orion: Aso et al. 2000). 

Since possible evolutionary processes of stable globules are expected not to be dynamical, their 
lifetime should be long compared with unstable ones. However, in our starless globule samples, 
the fraction of stable globules (density contrast is less than 14) is small (3 out of 11 sources). The 
observed starless globules are mostly populated near the critical Bonnor-Ebert sphere, and decreases 
toward the smaller density contrast (Figure 6). This result is inconsistent with an idea that all of 
nearly critical globules originated from stable low-density-contrast globules; if the idea is the case, 
the number of the observed low-density-contrast globules should be larger than or comparable with 
that of the nearly critical ones. Thus, it is most likely that most of the nearly critical globules 
which we observed are not evolved from stable globules. The origin of nearly critical globules can 
be related to the formation mechanism of globules. Further studies on this subject are desirable. 

5.5. Slow Collapse of A Nearly Critical Bonnor-Ebert Sphere 

After gravitational collapse begins, globules evolve toward higher central condensation. If 
we assume an extreme case in which all the globules fitted as unstable Bonnor-Ebert equilibrium 
states are already collapsing, the difference in £ max is attributed to different evolutionary stages of 
globules. Since the Bonnor-Ebert spheres provide a good fit to starless globules, it is meaningful to 
consider the gravitational collapse of the Bonnor-Ebert sphere (e.g., Foster & Chevalier 1993; Ogino, 
Tomisaka, & Nakamura 1999). We expect that the critical Bonnor-Ebert sphere can approximate 
the initial condition of collapse, because a large fraction of the starless globules are located near 
the critical state, £ max = 6.5 ± 2, as shown in § 5.1. 

In order to make quantitative comparisons of the density structures between our observations 
and a collapsing gas sphere starting from a nearly critical Bonnor-Ebert sphere, we used the result 
of theoretical calculations from Aikawa et al. (2005). In their calculation, the central density n c (H) 
and kinetic temperature T of the critical Bonnor-Ebert sphere were set as n c (H) = 2 x 10 4 cm~ 3 
and T = 10 K, respectively. No cloud rotation, magnetic fields, or non-thermal turbulence were 
considered, and a fixed boundary condition of zero velocity was adopted at the outermost radius of 
R = 0.2 pc, which is slightly larger than the critical radius of an equilibrium sphere (see, Aikawa et 
al. 2005). We note that the initial condition parameters of n c (H), T, and R, are consistent with the 
observed physical properties of globules around £ max = 6.5 (see, Table 4 and 5). The initial density 
distribution was multiplied by a constant factor a of 1.1 in order to initiate collapse. The slight 
density enhancement of a = 1.1 is preferable to approximate the collapse from a nearly critical 
Bonnor-Ebert sphere. The time scale of the collapse is several times the free-fall time, ~ 1.2 x 10 6 
yr (~ 4 x tfi), which is comparable with the observational lifetime of starless dense cores (~ 10 6 yr; 
e.g., Visser et al. 2002; Lee & Myers 1999). 

Figure 8 shows the evolution of a radial column density profile of the collapsing sphere for the 
case of a = 1.1 described above. The gray solid line denotes initial column density profile (t = 
yr), and the black solid, dotted, and broken lines denote density profile at specific times of 4.4 x 10 5 , 
8.9 x 10 5 , and 1.1 x 10 6 yr, respectively. We found that the column density profiles of the collapsing 
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sphere are well fitted with the static Bonnor-Ebert spheres (equilibrium solutions) in the range of 
Cmax < 25 within the relative difference of 2%. The plus symbols show the best-fitting result for 
each collapsing profile with the Bonnor-Ebert sphere, which are labeled with corresponding £ max 
values; for example, the density profile of the collapsing sphere at 4.4 x 10 5 yr mimics the unstable 
Bonnor-Ebert sphere of £ max = 7.6. This result suggests that observed starless globules with large 
density contrasts are gravitationally collapsing objects. Figure 9 shows the relationship between 
Cmax versus density contrast. The solid line is for the Bonnor-Ebert sphere, and the dotted line 
is for the collapsing sphere which plots the best-fitting £ max values against the density contrast of 
the collapsing sphere. The plus symbols denote the elapsed time of the collapsing sphere at 10% 
intervals of the total collapse time. Since the free-fall time scale increases with decreasing density, a 
collapsing gas sphere starting from a nearly critical Bonnor-Ebert sphere should mimic marginally 
unstable Bonnor-Ebert spheres for a long time. Thus, this is qualitatively consistent with the fact 
that the £ max for the starless globules is mostly populated near the critical state (7 out of 11 sources) 
and the remaining starless globules show clearly unstable states (£ m ax > 10; 4 out of 11 sources) 
as shown in § 5.1. From the model prediction, we note that 50% of collapsing globules should be 
located in £ max = 6.5 — 8.5, and the rest of the globules should have £ max larger than 8.5. The 
fraction is consistent with our result that half of the starless globules with £ max > 6.5 are located 
in £max = 6.5 — 8.5 (4 out of 8 sources, see Fig. 5). In Figure 10, we demonstrate the histogram 
of logarithmic density contrast for starless globules (the same as Figure 6) and the time evolution 
of the model collapsing sphere (solid line with dots). Each dot represents the residence time of a 
collapsing sphere measured at 0.2 intervals in the logarithmic density contrast. The histogram of 
the observed globule in the unstable region is consistent with that from model calculations of a 
collapsing sphere. 

As discussed above, the density structure of the observed globules can be explained in terms 
of the slow collapse of a nearly critical Bonnor-Ebert sphere. The anticipated inward velocity from 
the model calculation remains subsonic^ for a long time; the peak inward velocity is less than ~ 0.1 
km s _1 for ~ 10 6 yr, which is ~ 80% of the total collapse time (see, Fig. 1 of Aikawa et al. 2005). 
Previous reports on the observations of CS (J = 2 — ► 1) (blue) asymmetry profile, indicative of 
gravitational collapse in dense cores, show that infall candidates constitute 25 % of the starless core 
samples (Lee, Myers, & Tafalla 1999). It is possible that the rest of the starless cores have a slow 
collapse motion which is hard to detect. Radio molecular line observations of infall asymmetry 
profiles with sufficiently high velocity resolution (e.g. ~ 0.01 km s^ 1 ) are necessary to confirm 
the scenario of slowly collapsing globules. Among the starless globules in Figure 5, the gas inward 
velocity for Lynds 694-2 (£ max = 25 ±3) was measured as 0.05 — 0.07 km s _1 (Lee, Myers, Sz Tafalla 
2001; Lee, Myers & Plume 2004), which is consistent with the prediction from the collapse model 
(a = 1.1). 



^In § 4.2.2, we have estimated distance to the globules using T e g from radio molecular line observations as shown 
in Table 3. A collapsing motion in the globules should contribute to the derivation of T e g if the motion is large 
enough compared with the observed molecular line width. Since predicted inward velocity from the model is small, 
it is likely that the contribution of collapsing motion to the derived distance is negligible. 
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6. Summary 

We have carried out a near-infrared (J, H, and K s ) imaging survey of ten Bok globules using 
the infrared camera SIRIUS on the IRSF 1.4 m telescope in South Africa. Ay distributions for 
the ten globules were derived through measurements of the H — K color excess or the stellar 
density at J, H, and K s (star count). The radial column density profile for each globule was 
analyzed using the Bonnor-Ebert sphere model, which well fits the profile within observational 
uncertainties, and physical properties of globules, e.g., size, central density, temperature, mass, 
external pressure, and center-to-edge density contrast, were derived. The dimensionless radial 
parameter £ max determines the shape of a Bonnor-Ebert density profile as well as the stability of 
the equilibrium state, e.g., the solution of £ max > 6.5 is unstable to the gravitational collapse. We 
investigated the stability of globules on the basis of £ max for ten globules from our observations and 
four globules in the literature. In addition to the near-infrared imaging, we have carried out radio 
molecular line observations toward the same ten Bok globules using the 25-element focal-plane SIS 
receiver BEARS on the 45 m telescope of Nobeyama Radio Observatory (NRO). We measured line 
width of the globules for independent measurements of the "effective" temperature including the 
contribution of turbulence. It was confirmed that most of the globules are dominated by thermal 
support. Since the distance-dependent parameters of the Bonnor-Ebert sphere were derived using 
the initially assumed distance to the globules, derived best-fitting temperature was affected by the 
distance assumption. We constrained the distance to the globules by comparing the temperature 
from the molecular line width with the temperature from the Bonnor-Ebert fitting, and rescaled 
distance-dependent physical quantities of globules. 

(1) We found that more than half of the starless globules (7 out of 11 sources) are located 
near the critical state, £ max = 6.5 ± 2. Thus, we suggest that a nearly critical Bonnor-Ebert 
sphere characterizes the typical density structure of starless globules, and it approximates the 
initial condition of gravitational collapse. 

(2) We found that four out of eleven starless globules show clearly unstable states (£ max > 10). 
Since unstable equilibrium states should not be long sustained, we expect that they are already 
collapsing toward higher central condensation or that extra force (e.g., magnetic and/or turbulent 
pressure) accounting for large £ max stabilizes the globules. It was also found that all three star- 
forming globules have unstable solutions of £ max > 10, which is consistent with the fact that they 
have started gravitational collapse. 

(3) We investigated the collapse of the Bonnor-Ebert sphere from a nearly critical state using 
the model calculation of Aikawa et al. (2005), and found that the column density profiles of 
the collapsing sphere mimic those of static Bonnor-Ebert spheres (unstable equilibrium solutions). 
By relating £ max to the collapsing sphere at a specific time, the evolutionary state of globules 
can be interpreted, and the detection probability of each £ max value can be predicted from the 
model calculation. Since the evolutionary timescale decreases with increasing density, the collapsing 
sphere resembles a marginally unstable Bonnor-Ebert sphere for a long time. It was found that the 
frequency distribution of £ max for the observed starless globules is consistent with that from model 
calculations of the collapsing sphere. 
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Appendix A: Possible Systematic Uncertainties in £ max 

Here we discuss the ambiguity of derived column density (Ay), which leads to the systematic 
uncertainties in £ max - In order to evaluate the dependence of the conversion factor, Ay / Eh-k 
and Ay I A\ (A = J, H, and K), on the ratio of total-to-selective extinction Ry, we used the in- 
dependent extinction law empirically derived by Cardelli et al. (1989). Though it is well established 
that Ry has a uniform value of ~3.1 in the ISM (e.g., Whittet 1992), there are some observational 
reports on higher Ry values up to ~ 6 in dense molecular clouds (e.g., Kandori et al. 2003, and 
references therein). The empirical extinction law for Ry = 6 gives the conversion factors smaller 
~ 20 % than those for the case of Ry = 3.1, which we adopted in the analysis (§3.1 and §3.2). 
Note that the extinction law gives slightly different conversion factors compared with Rieke & 
Lebofsky (1985) as shown in Table 3 of Cardelli et al. (1989). This possible systematic uncertainty 
as well as the uncertainty of the gas-to-dust ratio can be included in the derivation of A(H2). We 
have converted Ay to H2 column density using the relationship iV(H2)/Ay = 9.4 x 10 20 cm™ 2 
mag™ 1 derived from Bohlin, Savage, Sz Drake (1978). Similar relationships were obtained from 
the X-ray absorption studies toward SNRs with high dynamic range (Ay < 30 mag) as Ah /Ay = 
1.8 - 2.2 x 10 21 cm™ 2 mag™ 1 (Gorenstein 1975; Ryter et al. 1975; Ryter 1996; Predehl k Schmitt 
1995). Thus, the conversion coefficient is likely to be determined within ~ 20 % error. It can be 
seen from the panel (a) of Figure 7 that Ay is roughly correlated with about l/2th power of density 
contrast for the Bonnor-Ebert spheres. Thus, for the case that Ay is systematically changed by 
a factor the value of density contrast varies by /3 L5 accordingly. Since £ max is proportional to 
y/fc, corresponding £ max value is expected to vary by /3°- 75 . If we assume that derived Ay for 
each globule is overestimated by 20 %, £ max is overestimated by ~ 15 %. This is comparable to 
the uncertainty in £ max derived from the \ 2 fitting, and this modification should not significantly 
change our conclusions on the density structure of globules. 
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Table 1: Source list 



Name 


R.A. (J2000) a Dec. (J2000) a 


Distance 


IRAS b 


Other name 


Reference 




m 


(° ' ") 


(PC) 








CB 87 


17 25 05 


-24 07 19 


160 


No 


Barnard 74, Lynds 81 


1,2 


CB 110 


18 05 55 


-18 25 10 


180 


No 


Lynds 307 


1,3,4 


CB 131 


19 17 00 


-18 01 52 


180 


No 


Barnard 93, Lynds 328 


1,3,4 


CB 134 


18 22 45 


-01 42 40 


260 


No 




1,5,A 


CB 161 


18 53 56 


-07 26 29 


400 


No 


Barnard 118, Lynds 509 


1,6 


CB 184 


19 31 52 


+16 27 14 


300 


No d 


Lynds 709 


1,3,4 


CB 188 


19 20 16 


+11 36 15 


300 


Yes e 


Lynds 673-1 


1,3,4 


FeSt 1-457 


17 35 45 


-25 33 11 


160 


No 




1,3,4,B 


Lynds 495 


18 38 58 


-06 44 00 


200 


No 




1,7,C 


Lynds 498 


18 40 11 


-06 40 45 


200 


No 




1,7,C 



"Position of field center for each near- infrared image. 

^Existence of IRAS point sources within optical boundary of each core. 

c (l) Dutra k Bica (2002); (2) Huard, Sandell & Weintraub (1999); (3) Launhardt & Henning (1997); (4) Dame et al. (1987); (5) Straizys, Cernis, & Bartasiute 
(1996); (6) Leung, Kutner, & Mead (1982); (7) Schneider, & Elmegreen (1979); (A) Assumed to have the same distance as Serpens molecular cloud; (B) Assumed 
to have the same distance as Barnard 83; (C) Assumed to have the same distance as GF 5 dark cloud filament. 

d IRAS 19116+1623 (Class II located near core boundary): Launhardt (1996) 

e IRAS 19179+1129 (Class I with outflow): Launhardt (1996), Yun, & Clemens (1992) 



Table 2: Results of stellar detection 



Field Limiting mag. a Threshold mag. b Detected number of stars c 



J H K s J H K s J H K s H and K, 



CB 87 


18. 


,9 


18.4 


17.5 


18.5 


17.8 


17.1 


8107 


10588 


6510 


7609 


CB 110 


18. 


.8 


17.9 


17.1 


18.1 


17.2 


16.4 


12471 


14220 


12340 


12119 


CB 131 


19. 


1 


18.6 


17.5 


18.1 


17.4 


16.8 


6678 


10760 


111868 


9890 


CB 134 


19. 


.5 


19.2 


17.9 


18.9 


18.5 


17.5 


5631 


7530 


4530 


5510 


CB 161 


18. 


.5 


18.3 


17.3 


18.0 


17.7 


16.9 


6431 


8406 


5330 


6044 


CB 184 


18. 


.8 


18.1 


17.3 


18.3 


17.8 


17.1 


5235 


6056 


4241 


4673 


CB 188 


18. 


,9 


18.3 


17.4 


18.4 


17.6 


17.0 


6743 


8963 


7372 


7669 


FeSt 1-457 


19. 





18.2 


17.4 


18.6 


17.2 


16.8 


9261 


10614 


10444 


9240 


Lynds 495 


19. 


.1 


18.0 


16.6 


18.4 


17.2 


15.9 


10764 


13404 


10510 


10760 


Lynds 498 


18. 


.7 


17.6 


16.7 


17.7 


16.8 


16.0 


13826 


13696 


12064 


11673 



"Limiting magnitude is defined by the magnitudes whose photometric error is equal to 0.1 mag. 

fc Threshold magnitude is defined by the peak value of the histogram distribution of magnitudes for detected stars whose magnitude errors are less than 0.1 mag. 
c The first three columns show the number of stars detected on the J, H, or K s images (single band detection) whose photometric error is less than threshold 
magnitudes. The last column shows the number of stars detected both on the H and K 3 images whose photometric errors, sj SH' 2 + SKg, are less than 0.15 mag. 



Table 3: Molecular line parameters of globules 



Name 


Vlsr 


AV a 


AF NT b 


T c g c 


r d 


rp C 


Molecular Line 




(km s _1 ) 


(km s _1 ) 


(km s _1 ) 


(K) 




(K) 








CB 87 f 


4.73±0.02 


0.21±0.04 


0.17±0.05 


11.4±0.8 


0.55±0.12 


10 


C 18 (J = 


1 


-0) 




4.73±0.01 


0.25±0.02 


0.22±0.03 


12.4±0.6 


0.38±0.04 


10 








CB 110 


5.94±0.02 


0.50±0.04 


0.48±0.04 


21.8±2.2 


0.62±0.06 


10 


C 18 (J = 


1 


-0) 




5.91±0.01 


0.45±0.02 


0.43±0.02 


19.3±1.0 


0.47±0.02 


10 








CB 131 


6.66±0.01 


0.56±0.02 


0.55±0.02 


25.1+1.3 


0.58±0.03 


10 


C 18 (J = 


1 


-0) 




6.65±0.01 


0.57±0.02 


0.56±0.02 


25.8±1.1 


0.56±0.02 


10 








CB 134 


8.56±0.01 


0.28±0.02 


0.25±0.02 


13.2±0.6 


0.94±0.10 


10 


C 18 (J = 


1 


-0) 




8.562±0.004 


0.24±0.01 


0.20±0.01 


12.1±0.2 


0.81±0.04 


10 








CB 161 


12.48±0.01 


0.25±0.02 


0.22±0.02 


12.5±0.6 


0.61±0.06 


10 


C 18 (J = 


1 


-0) 




12.47±0.01 


0.27±0.01 


0.24±0.02 


13.0±0.4 


0.44±0.02 


10 








CB 184 


6.10±0.01 


0.35±0.02 


0.33±0.02 


15.5±0.6 


1.01±0.08 


10 


C 18 (J = 


1 


- 0) 




6.096±0.005 


0.35±0.01 


0.33±0.01 


15.5±0.4 


0.88±0.04 


10 








CB 188 


6.98±0.01 


0.44±0.03 


0.42±0.03 


19.0±1.3 


0.54±0.04 


10 




i 
± 


— m 




a Q9-i-n m 






oo q_i_i 9 




1 n 








FeSt 1-457 


5.820±0.003 


0.182±0.006 


0.132±0.009 


10.9±0.1 


9.56±1.29 


6.25±0.21 


N 2 H+ (J = 


1 


-0) 




5.824±0.002 


0.175±0.004 


0.122±0.005 


10.8±0.1 


11.48±0.86 


5.80±0.09 








Lynds 495 


12.18±0.10 


0.26±0.02 


0.22±0.02 


12.6±0.6 


0.60±0.06 


10 


C 18 (J = 


1 


-0) 




12.18±0.01 


0.26±0.01 


0.23±0.02 


12.7±0.4 


0.54±0.03 


10 








Lynds 498 


12.51±0.01 


0.18±0.03 


0.14±0.04 


11.0±0.6 


0.41±0.07 


10 


C 18 (J = 


1 


-0) 




12.48±0.01 


0.29±0.02 


0.27±0.03 


13.6±0.7 


0.31±0.02 


10 









a Observed FWHM line width. 

b Non-thermal line width. We assumed thermal line width AVt is equal to 0.443 km s _1 (Tk = 10 K). 
c Temperature including non-thermal (turbulent) contribution (T e ft = Tk + Tnt). 

d Peak optical depth. The value for FeSt 1-457 is the sum of the peak optical depths of the seven hyperfine components. 
e We assumed the excitation temperature for the C O spectra to be 10 K. 

^The first row refers to the values from the spectrum toward the center of each globule, and the second row refers to the values from the composite spectrum of the 
central nine positions. The error represents la error in the line profile fitting. 



Table 4: Physical properties of globules I 



Name 


R.A. (J2000) a Dec. (J2000) a 


Br 


£max 


^c/ ^edge 


A v c 


e d 


P.A. d 


Tfit 6 


x 2 


Stability f 


Reference g 




/h m 


/ ft \ 

(° ) 


(") 






(mag) 




(°) 


(") 








CB 87 


17 24 58.1 


-24 06 49 


87.5+2.6 


5.1+0.4 


8.1 


6.6 


0.52 


120 


<110 


3.19 


Stable 


A 


CB 110 


18 05 54.7 


-18 25 02 


61.1+3.1 


14.0+3.0 


94.3 


19.7 






<70 


1.17 


Unstable 


A 


CB 131 


18 17 00.5 


-18 02 05 


103+7.4 


16.3+5.1 


139 


32.9 








<150 


3.22 


Unstable 


A 


CB 134 


18 22 46.3 


-01 42 50 


59.6+4.3 


18.5+4.9 


187 


26.8 






<55 


0.44 


Unstable 


A 


CB 161 


18 53 56.0 


-07 26 07 


62.5+6.9 


8.1+1.4 


25.1 


13.0 


0.72 


125 


<50 


0.11 


Unstable 


A 


CB 184 


19 13 51.7 


+16 27 27 


112+20 


8.1+1.6 


24.9 


11.2 


0.68 





<75 


0.15 


Unstable 


A 


CB 188 


19 20 15.7 


+11 36 07 


127+19 


16.0+2.9 


132 


27.2 






<100 


5.17 


T T j_ 1 1 

Unstable 


A 


FeSt 1-457 


17 35 47.5 


-25 32 59 


144+11 


12.6+2.0 


74.5 


41.0 






<140 


1.41 


Unstable 


A 


Lynds 495 


18 38 57.4 


-06 44 06 


75.0+7.8 


7.2+1.4 


18.6 


12.7 






<65 


0.14 


Unstable 


A 


Lynds 498 


18 40 10.5 


-06 40 43 


75.0+3.0 


4.7+0.4 


6.63 


10.9 






<70 


0.84 


Stable 


A 


Barnard 68 






100 


6.9+0.2 


16.6 












Unstable 


B 


Barnard 335 


19 37 00.9 


+07 34 10.0 


125 


12.5+2.6 


73.1 








<100 


2.93 


Unstable 


C 


Coalsack 


12 31 38.6 


-63 43 42.5 


290 


5.8 


10.9 


11.5 






55-400 




Stable 


D 








140 


7.0+0.3 


17.2 


23.0 










Unstable 


E 


Lynds 694-2 


19 41 04.4 


+10 57 00.9 


54 


25+3 


364 








<83 


1.12 


Unstable 


F 



"Position of core center (centroid of Ay distribution). 

6 Center-to-edge density contrast determined from £ ma x value. 

c Ay at core center. 

d Fitted ellipticity and position angle. 

e Profile fitting region in the Bonnor-Ebert fit. 

■^Stability of equilibrium state of the Bonnor-Ebert sphere. 

9 (A) This work; (B) Alves et al. (2001); (C) Harvey et al. (2001); (D) Lada et al. (2004); (E) Racca et al. (2002); (F) Harvey et al. (2003) 
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Table 5: Physical properties of globules II 







D h 


R c 




77 C 




M c 


p . 

1 cxt 


c 


R ofovovi CO 

x l A~±r~i c-J-i^v; 




(K) 


(vc) 


(AYS) 




(cm - 


I 




(K cm 


I 




CB87 


6.0+0.4 


160 


(1 40+0 04' 


xlO 4 


(3.2+0.4' 


xlO 4 


76+0 07 


(2.9+0.4' 


xlO 4 


A 




1 1 4 

X X ."-i 


304 

Out 


(n Rfi+o 08' 


xlO 4 

/\ XU 


fi 7_i_o 9' 


xlO 4 


o 7Q+0 94 


in q+o 4^ 


xlO 4 

A 1U 


A 


CB 110 


7 n_i_n a 


1 80 

XOU 


(-t i o+O 


xlO 4 


(a 5+1 7' 


x 10 5 


74+0 1 7 


(A o+O 


xlO 4 


A 




21.8 


561 


(O AQ-Lfi 1B> 

\ *J .'rtJ—l—L/. X \J j 


xlO 4 


fi 5+0 fi A 


x 10 5 


7.21+1.64 


(A o+O 9' 


xlO 4 


A 


CB 131 


14.0±0.8 


180 


(I 85+0 13' 


xlO 4 


(4.4+2.3) 


xlO 5 


2 43+0 73 


(5.3+1.5) 


xlO 4 


A 




25.1 


323 


(3 32+0 24' 


xlO 4 


(2.5+1.3' 


xlO 5 


7 83+2 35 


C5 o+i c-| 


xlO 4 


A 


CB 134 


1 Q niri o 

XO . \J—1—\J . o 


260 


\ X . ' " l I Uill i 


xlO 4 


C7 r+q 4 N 


x 10 5 
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"First row: value from the Bonnor-Ebert fit using assumed distance to the globules (see Table 1), Second row: value 
from molecular line observations (T e tf; see, Table 3). 

b First row: same value as listed in Table 1, Second row: calibrated value by using T e g to scale the fitted Bonnor-Ebert 
model parameters (see, §4.2.2). 

c The values in the first and second rows are derived using the distances listed in the third column. 
d (A) This work; (B)-(F) Same references as those in Table 4; (G) Hotzel et al. (2002a). 




Fig. 1. — J, H , and K s three-color composite images for ten Bok globules (J: blue, H: green, K s 
red). Each image size is 930x930 pixel (~7 arcminutes). North is up, and east is toward the left. 
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Fig. 2. — Ay distribution for Bok globules. The area enclosed by the broken line is the reference 
field expected to be free from dust extinction, and the area enclosed by the dotted line is the masked 
region which is excluded in the derivation of the radial Ay profile. 
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Fig. 2. — continued. 
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Fig. 3. — Line profiles toward the center of ten globules (N2H + for FeSt 1-457, C 18 for the other 
globules). The gray solid line denotes fitting result of each spectrum. 
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Fig. 4. — Radial column density profile and the Bonnor-Ebert model fit for the ten Bok globules. 
The dots and error bars represent the average iV(Ha) values at each annulus arranged at intervals 
of 9" and the rms dispersion of data points in each annulus, respectively. The solid line denotes 
the best-fit Bonnor-Ebert profile which is convolved to the observing resolution (typically ~30") 
with the same beam as used in the Ay measurements. The dashed line denotes the Bonnor-Ebert 
profile before the convolution. The gray dot-dashed line denotes the 1 a deviation of the column 
density measured in the reference field. 



- 34 - 



CB 106 F&5! 1-45? 




Lynda 495 Lynds 4S& 




Fig. 4. — continued. 
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Fig. 5. — The solid line shows the relation between £ max an d density contrast (center-to-edge 
density ratio). The location of our 10 globules and previously reported 4 globules (Barnard 335, 
Harvey et al. 2001; Barnard 68, Alves et al. 2001; Coalsack Globule II, Racca et al. 2002; Lada et 
al. 2004; Lynds 694-2, Harvey et al. 2003) are shown in the figure. Filled and open black circles 
denote starless and star-forming globules, respectively. Filled gray circles denote the fitting results 
when the masked region on the Ay map (see Fig. 2) is included in the derivation of the column 
density profile. The vertical dashed line denotes critical state of £ m ax= 6.5. In the lower panel of 
the figure, horizontal error bars denote la error in the x 2 density profile fitting. 



-36- 



] Starless [ ] Star -forming 



5 

2 
1 


5 



i: 



-n — i 1 1 r 

Critical 



0.5 1 1.5 2 2.5 3 3.5 

Log (Density Contrast) 



Fig. 6. — Histograms of logarithmic density contrast for starless and star-forming globules. The 
vertical broken line denotes the density contrast of the critical Bonnor-Ebert sphere. 
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Fig. 7. — Correlations between density contrast and the other physical parameters. The panels of 
(a) to (f) show the correlation diagrams for Ay, central density, radius, mass, temperature, and 
external pressure, respectively. Filled and open circles denote starless and star-forming globules, 
respectively. The vertical gray line denotes the density contrast of the critical Bonnor-Ebert sphere. 
The plotted data points are the physical quantities after distance-correction (see, §4.2.2, Table 5). 
The physical quantities without distance-correction are shown as the symbol "— " in each panel. The 
broken lines represent the relationships for the Bonnor-Ebert spheres with a constant temperature 
T cff of 13.2 K and P cxt of 5.3 x 10 4 K cm" 3 . 
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Fig. 8. — Radial column density distribution of a collapsing gas sphere at some specific times (black 
solid, dotted, and broken line) started from the initial condition of a nearly critical Bonnor-Ebert 
density structure (gray solid line). Dotted and broken lines are shifted by 5 and 10 (10 21 cm -2 ) in 
the figure, respectively. Best-fit solution of the Bonnor-Ebert sphere for each line is also plotted as 
plus symbol. 
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Fig. 9. — The solid line shows the relation between ^ max and density contrast for the Bonnor-Ebert 
sphere (same as Fig. 5). The dotted line denotes the similar relation for a collapsing sphere starting 
from a nearly critical Bonnor-Ebert sphere (a = 1.1 model; Aikawa et al. 2004); the column density 
profiles of the collapsing sphere were fitted with static Bonnor-Ebert sphere solutions, and the best- 
fit £max values are plotted against the density contrast values of the collapsing sphere. The plus 
symbols denote the elapsed times from the onset of collapse at 10% intervals of the total collapse 
time (t co i). 
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Fig. 10. — The histogram of the logarithmic density contrast for starless globules (the same as 
Figure 6). The solid line with dots shows the time evolution of the model collapsing sphere; each 
dot represents the residence time of a collapsing sphere measured at 0.2 intervals in the logarithmic 
density contrast. The vertical axis on the right side shows the percentage of the residence time 
with respect to the total collapse time. 



